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Abstract 

We investigated the effect of time delays on phase configurations in a set of 
two-dimensional coupled phase oscillators. Each oscillator is allowed to inter- 
act with its neighbors located within a finite radius, which serves as a control 
parameter in this study. It is found that distance-dependent time-delays in- 
duce various patterns including traveling rolls, square- like and rhombus-like 
patterns, spirals, and targets. We analyzed the stability boundaries of the 
emerging patterns and briefly pointed out the possible empirical implications 
of such time-delayed patterns. 
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Spatiotemporal patterns arise in numerous physical, chemical, and biological systems [1]. 
The brain, one of the most complex systems, is now also known to generate spatiotemporal 
patterns such as plane waves and spirals [2-6]. Early studies, in this context, have concen- 
trated on exploring the dynamics arising in a set of coupled oscillators [6-11]. It has been, 
however, generally assumed that interactions among the individual oscillators are instanta- 
neous. Recently, noting that an inclusion of time delay is more natural in realistic systems, 
several authors have investigated the effects of time delay, and have found the consequent 
dynamic phenomena such as multistability, desynchronization, clustering, amplitude death, 
anticipated synchronization and slow switching [9-16]. It has been also found recently that 
a time delay, proportional to a distance between elements, could produce a propagating 
structure in a one- dimensional coupled-oscillators [17]. 

In this paper, we discuss further the effect of time delays on the dynamics in a complex 
system in two dimensions, employing a two-dimensional coupled oscillator model [2-6]. It is 
found that time delays among constitutive elements alone can induce definite spatiotemporal 
patterns in two-dimensional complex system. 

We start with the following coupled oscillator model equations for phase variables 6^ (t) 
(i and j are integers) : 



jy- 0<r k i tij <r o 
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0<r k i }i j<ro 

with N(r ) = W(r k nj). Here r k nj denotes the distance between element 

k,l 

and element (k,l). The metric properties of the ensemble are fixed by setting r^jj = 
- k) 2 + (j - I) 2 . In this model an interaction is non-local and time delayed, character- 
ized by coupling length r and signal propagation speed v, and time delays are expressed 
by fki t ij/v. Assuming that an individual element interacts with its neighbors located within 
a finite radius, we introduce a coupling radius denoted by r ; We also discussed the case of 
r — > oo in this study. Further we introduce a weighting function W{rki,ij) which mimics 
the interaction strength depending on a distance between interacting elements. Since there 
is no law about W(rki,ij), as far as we know, we choose W(rki,ij) = l/rki,ij assuming that 



the Gauss' law of interactions in two dimensional cases is also valid here [18]. 

The above equation without time delays is a model of coupled limit cycle oscillators in 
the so-called phase approximation [7]. It is known that, in the absence of delays rkijj/v , 
for any positive values of the coupling coefficients K, all the oscillators finally form a planar 
solution where all the oscillators are synchronized, oscillating with frequency equal to that 
of the individual oscillators [7]. The present study of the model (1), on the other hand, 
shows that the final steady states correspond to planar solutions with frequencies different 
from w or to frequency-synchronized nontrivial phase configurations. 

To find out possible phase configurations, let us assume steady solutions in which all the 
oscillators have the same frequency f2. Then the solutions of the model can be written as 
9ij(t) —Qt + 4>ij. Substitution of these solutions into Eq. (1) yields 

j£ 0<r fc!j4j <r j 

fi = w + ati \ Yl sin(0 H - <f>ij - Vtr khij /v). (2) 
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0<rki,ij<ro 

This self-consistency condition requires that the sums N ^ sin(0 fc / — 0^ — 

^ r ki,ij/ v ) be same for all [17]. Let us consider the situation in the continuum limit. 

Differentiating the term with respect to f\— we get the following. 

^0 = ^?, for all rand P. (3) 
dr' df w 

This implies that, in the continuum limit, <f>(f) = 7r • e + O is a unique solution, up to 

the arbitrariness of O an d e, as long as e is a unit vector embedded in the lattice plane. 

The case of 7 = corresponds to synchronized planar solutions, while nonzero 7's make roll 

solutions traveling in the direction of e with the phase velocity 

However, the space discreteness of the system admits topological defects, so other phase 

configurations are also possible [19]. In this study, we set w = ir/5, and solved the equation 

on a 64 by 64 square plane (1 < i < 64, 1 < j < 64) with periodic boundary conditions. We 

checked that different size of lattice does not alter the result. We tried several kinds of initial 

conditions. As the initial conditions in phase of the form ^-,(0) = 9 + in this study, 

where 77^ denotes a random number, we considered two different cases: (i) < rjij < 2n, 



and (ii) < rjij < 10~ 6 . In addition, for the dynamic history of individual oscillators for 
t < 0, we considered two situations (i) where the oscillators were stationary, and (ii) where 
the oscillators have evolved independently from one another with their natural frequency w. 

Figure 1, which is made with the parameters K = 0.6 and r = 10, shows effects of 
the delay l/v. When the delay l/v is small (< 0.6), only planar solutions^ = 0) occur. 
When the delay is large enough(> 1.1), planar solution is not available but patterns with 
nonzero 7's occur. In the intermediate region (0.6 < l/v < 1.1), planar solutions as well 
as patterns are possible. The lines in Fig.l are the numerically obtained solutions of Eq. 
(2) with (fiij = ji, 7 = 0, ||, ||, and ||. We see that the wavelengths of patterns become 
shorter generally as delay increases. 

Figure 2 provides the phase diagrams of the model. If the coupling strength K is above 
K c2 (denoted by open squares), only a synchronized planar solution appears, but if K is be- 
low the critical value K c \ (denoted by filled squares), planar solution is not available and the 
system evolves into patterns such as rolls [as in Fig. 3(a)], or other phase distributions [as 
in Figs. 3(b) - 3(d)] depending on initial configurations. In the regions of K c \ < K < K c2 , 
we see both patterns and planar solutions depending on initial configurations. Figure 4 
shows the patterns that appear only in the limited regions in the phase diagram. When 
r is relatively small (^ 5), in other words, when the effect of space discreteness becomes 
important, spirals occur more frequently [See Fig. 4(a)] than other configurations. Further- 
more, the model sometimes organizes target patterns from slightly perturbed homogeneous 
phase configurations for parameter K ~ K c \ [see Fig. 4(b)]. With K ~ K c2 , there emerge 
patterns embedded in a planar solution as shown in Fig. 4(c). 

Note that the arrangement of the phase singularities plays an essential role in forming 
the patterns of Figs. 3(b) - (d). The singularities are possible due to the discreteness of 
the space as in spiral patterns, but are dynamically arranged in a rhombus [Fig. 3(b)] 
or a square lattice [Fig. 3(c)]; the phases of oscillators make rotating waves around those 
singularities [20]. However, rolls are more prevalent than squares or rhombus patterns in 
the entire region of the diagram. 
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There seem to exist three kinds of multistabilities in our model. First, as stated above, 
we can see planar solutions, rolls or other phase distributions for the same parameters 
depending on initial configurations [see Figs. 1-3]. Second, there are multiple states with 
different wavelengths and frequencies even for the same parameters in some cases, which is 
why we see multiple Si's of patterns in Fig. 1. This kind of multistability is also reported in 
one-dimensional lattices [17]. Third, different initial configurations make the different angles 
between two primitive translation vectors of patterns like Figs. 3(b) and 3(c). It seems that 
multistabilities are the major characteristics of time delayed systems [12-14]. 

In the systems of two-coupled oscillators [9,12], or globally coupled oscillators with uni- 
form time delays [13], the dynamics are determined by the value of fir, where Q is the 
frequency of the synchronized state and r is a time delay between oscillators. The physical 
meaning of the value is clear. It denotes the virtual phase difference between oscillators due 
to delay. If it is near 2nir, where n is an integer, each oscillator can keep pace with others to 
make the whole system synchronized, otherwise it cannot [13]. Likewise, we suggest that the 
dynamics of our model is determined by the quantity Qr /v, which we call 6 here, where Q 
is a frequency of a synchronized planar solution. Our conjecture is that there exists some 
value O c so that planar solutions are possible only for < 9 C [21]. 

We measured C directly from numerical simulations using Eqs. (1) and (2) in Fig. 5. 
We see that 6 C is a nearly constant value (~ 2.26), from which we can calculate K cl in 
Fig. 2. Substituting the roll solution = to the continuum limit of Eq.(2) makes the 
following equation. 

£1 = w / sin(rf2/w) J (7r)o?r, (4) 

r Jo 

where J is a Bessel function. In the case of planar solutions (7 = 0), the above equation 
can be integrated to give 

2sin 2 f V wr l J V ; 

The solid lines in Figs. 2 (a) and (b) are obtained by substituting = 2.26 to Eq.(5), and 
show good agreements with the simulation results. 



It is useful to look at what the existence of constant C implies. It means the following. 
Coupled oscillators tend to be synchronized, but when too large delays make the 6 value go 
beyond the tolerance limit 9 C , they fail to become a synchronized planar solution and take 
other configurations. The occurrence of incoherence in the coupled oscillators with uniform 
time delay [13,14] can now be explained in the same manner. 

It is also noteworthy that, in the limit of r — > oo, Eq. (5) converts to K c = jJ^rW, which 
is independent of 1/v. This is surprising, for it means that we must consider the effects of 
time delay however small when the interaction length is infinite as in most physical systems. 

To summarize, we have investigated the effects of time delayed interactions in an ensemble 
of coupled oscillators in two dimensions, and found that distance-dependent time delay 
alone can induce various spatial patterns, while without delay or with uniform delay alone 
only synchronized planar solutions might appear. We have analyzed also that the stability 
of planar solutions can be determined by the quantity Qr /v. Our study may find its 
relevance in recent experimental observations that non-locally interacting oscillators, which 
are especially probable for neurobiological systems, exhibit traveling rolls and spirals [4,6]. 
Finally, our results seem to suggest that time delay may play a significant role in the study 
of memory storage and information processing based on the spatiotemporal activities of 
neurons [2-6]. 
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Figure Captions 

Fig. 1. Synchronization frequency Q as a function of a delay (K = 0.6, r = 10); squares 
indicate planar solutions; circles indicate patterns. Lines are the solutions of Eq.(2). 
See the text for details. 

Fig. 2. The phase diagrams of the model with r = 10 in (a), and with l/v — 1 in (b). 
Above open squares, there is no pattern available, and only planar solutions appear. 
Between open and filled squares, the system evolves into planar solutions or patterns 
according to initial conditions. Below the filled squares, planar solution is not possible 
here, and only patterns appear. Solid lines are from Eq.(5). See the text for details. 

Fig. 3. Typical patterns generated in the model : (a) a traveling roll, (b) a rhombus-like 
pattern, (c) a square-like pattern, (d) coexistence of a rhombus-like pattern and a roll. 
K = 1.25, r = 25, and l/v — 1, for (a), (b), and (c). Note that planar solutions are 
also found for this set of parameter values.; K = 0.6, r = 10, and l/v = 1.1, for (d). 
See Fig. 4 (d) for the color coding. 

Fig. 4. Patterns obtained in simulations of the model : (a) spirals (K = 1.0, r = 4, and 
l/v — 1), (b) a target pattern (K = 0.4, r = 7, and l/v — 1), (c) a spiral embedded in 
a planar solution (K = 0.8, r = 8, and l/v = 1). In (c), the planar solution oscillates 
much more slowly than the spiral. Figure (d) shows the color code used in this paper. 

Fig. 5. The quantity Qr /v determines the stability of planar solutions. We cannot see 
planar solutions above the filled squares. 
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